#ifndef __MATRIX3_H__
#define __MATRIX3_H__

#include "Matrix.h"
#include "Vector3.h"

template <class Real>
class Matrix3 : public Matrix<3,Real>
{
public:
	// construction
	Matrix3();
	Matrix3( const Real entry[3][3] );
	Matrix3( const Matrix<3,Real>& matrix );
	Matrix3( Real f0, Real f1, Real f2, 
			Real f3, Real f4, Real f5,
			Real f6, Real f7, Real f8);
	
	Matrix3 ( const Real enteries[9], bool rowmajor );
	Matrix3 ( const Vector3<Real>& u, const Vector3<Real>& v,
		const Vector3<Real>& w, bool columns );
	Matrix3 ( const Vector3<Real>* v, bool columns );
	Matrix3 ( const Vector3<Real>& u, const Vector3<Real>& v );
	void MakeTensorProduct( const Vector3<Real>& u, const Vector3<Real>& v );

	// create a diagonal matrix
	Matrix3 ( Real m00, Real m11, Real m22 );
	void MakeDiagonal( Real m00, Real m11, Real m22 );

	// create rotation matrices ( positive angle - counterclockwise ),the
	// angle must be in radians, not degrees.
	Matrix3 ( const Vector3<Real>& axis, Real angle );
	void FromAxisAngle( const Vector3<Real>& axis, Real angle );

	// assignment
	Matrix3& operator= ( const Matrix3& m );
	Matrix3& operator= ( const Matrix<3,Real>& m );




	/*
	// operators
	Real* operator[] ( int row ) const;
	operator Real* ();
	Vector3 get_column( int col ) const;

	Matrix3& operator= ( int row ) const;
	bool operator== ( const Matrix3& matrix ) const;
	bool operator!= ( const Matrix3& matrix ) const;

	Matrix3 operator+ ( const Matrix3& matrix ) const;
	Matrix3 operator- ( const Matrix3& matrix ) const;
	Matrix3 operator* ( const Matrix3& matrix ) const;
	Matrix3 operator- () const;

	// matrix * vector
	Vector3 operator* ( const Vector3 vec ) const;

	// vector * matrix 
	friend Vector3 operator* ( const Vector3& vector,
        constMatrix3& matrix);

	// matrix * scalar
    Matrix3 operator* (Real fscalar) const;

	// scalar * matrix
    friend Matrix3 operator* (Real fscalar, const Matrix3& matrix);
*/

	// utilities
    Matrix3 Inverse() const;
	Matrix3 Adjoint() const;
	Real Determinant() const;

	// The matrix must be a rotation for these functions to be valid.
	void ToAxisAngle( Vector3<Real>& axis, Real& angle ) const;
	void Orthonormalize( );

	// The matrix must be orthonormal. The decomposition is yaw*pitch*roll.
	bool ToEulerAnglesXYZ( Real& yangle, Real& pangle, Real& rangle ) const;
	bool ToEulerAnglesXZY( Real& yangle, Real& pangle, Real& rangle ) const;
	bool ToEulerAnglesYXZ( Real& yangle, Real& pangle, Real& rangle ) const;
	bool ToEulerAnglesYZX( Real& yangle, Real& pangle, Real& rangle ) const;
	bool ToEulerAnglesZXY( Real& yangle, Real& pangle, Real& rangle ) const;
	bool ToEulerAnglesZYX( Real& yangle, Real& pangle, Real& rangle ) const;

	void FromEulerAnglesXYZ( Real yangle, Real pangle, Real rangle );
	void FromEulerAnglesXZY( Real yangle, Real pangle, Real rangle );
	void FromEulerAnglesYXZ( Real yangle, Real pangle, Real rangle );
	void FromEulerAnglesYZX( Real yangle, Real pangle, Real rangle );
	void FromEulerAnglesZXY( Real yangle, Real pangle, Real rangle );
	void FromEulerAnglesZYX( Real yangle, Real pangle, Real rangle );
	
	// SLERP (spherical linear interpolation) without quaternions.  Computes
    // R(t) = R0*(Transpose(R0)*R1)^t.  If Q is a rotation matrix with
    // unit-length axis U and angle A, then Q^t is a rotation matrix with
    // unit-length axis U and rotation angle t*A.
    static Matrix3 Slerp ( Real t, const Matrix3& r0, const Matrix3& r1 );

	void EigenDecomposition( Matrix3& rot, Matrix3& diag ) const;
	void SingularValueDecomposition( Matrix3& l, Matrix3& s, const Matrix3& r ) const;
	void SingularValueComposition( const Matrix3& l, const Matrix3& s, const Matrix3& r ) const;

	// factor M = Q*D*U with orthogonal Q, diagonal D, and upper triangular U
	void QDUDecomposition( Matrix3& q, Matrix3& d, Matrix3& u ) const;

	// special matrices
	static const Matrix3 ZERO;
	static const Matrix3 IDENTITY;

protected:
	
	// support for singular value decomposition
    static const Real epsilon;
    static const int max_iterations;

	// support for eigendecomposition
	void Tridiagonalize( Real diag[3], Real subdiag[3] );
	bool QLAlgorithm( Real diag[3], Real subdiag[3] );
    static void Bidiagonalize( Matrix3& ka, Matrix3& kl,
        Matrix3& kr );
    static void GolubKahanStep( Matrix3& ka, Matrix3& kl,
        Matrix3& kr );

    // support for spectral norm
    static Real max_cubic_root( Real coeff[3] );
};

typedef Matrix3<float> Matrix3f;
typedef Matrix3<double> Matrix3d;

#endif

